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A pair of oblique waves at low amplitudes is introduced in a supersonic flat-plate bound- 
ary layer. Their downstream development and the concomitant process of laminar to 
turbulent transition is then investigated numerically using Direct Numerical Simulations 
(DNS) and Parabolized Stability Equations (PSE). This abstract is the last part of an ex- 
tensive study of the complete transition process initiated by oblique breakdown at Mach 
3. In contrast to the previous simulations, the symmetry condition in the spanwise direc- 
tion is removed for the simulation presented in this abstract. By removing the symmetry 
condition, we are able to confirm that the flow is indeed symmetric over the entire compu- 
tational domain. Asymmetric modes grow in the streamwise direction but reach only small 
amplitude values at the outflow. Furthermore, this abstract discusses new time-averaged 
data from our previous simulation CASE 3 and compares PSE data obtained from NASA’s 
LASTRAC code to DNS results. 


I. Introduction 

To date, the most dominant nonlinear mechanism that eventually transitions a laminar, supersonic bound- 
ary layer to turbulence is still unknown. Knowledge of the relevant nonlinear mechanisms is however manda- 
tory for the accurate determination of transition onset. Previous investigations 1-3 of the nonlinear transition 
regime discovered two main nonlinear mechanisms, “oblique breakdown” and “asymmetric subharmonic res- 
onance”. Recently, a series of numerical studies 4-9 using direct numerical simulations (DNS) focused on 
several unresolved issues related to both mechanisms. In the first part of these studies, 4,5 the authors were 
able to identify oblique breakdown in the experiments by Kosinov and his co- workers, 2, 10-12 who investigated 
asymmetric subharmonic resonance initiated by a wave train in a flat-plate boundary layer at Mach 2. The 
identification of oblique breakdown in these experiments is of great importance since a detailed, experimental 
study of this mechanism has not been performed yet. 

The second part of these studies 8 addressed the question whether oblique breakdown or asymmetric 
subharmonic resonance, is a stronger nonlinear breakdown mechanism in supersonic boundary layers. In 
order to answer this question, Mayer et al. s investigated the early nonlinear transition regime initiated by 
a broad disturbance spectrum on a cone at Mach 3.5. To excite a wide range of disturbance waves, a 
wave packet was generated by a pulse through a hole on the cone surface. The disturbance spectrum in 
the frequency-azimuthal mode number plane of the wave packet exhibited traces of oblique breakdown and 
new features that could be explained with new resonance triads that are composed of three boundary layer 
modes with three different disturbance frequencies. Furthermore, it was shown that asymmetric subharmonic 
resonance is a limiting case of these new resonance triads with both secondary waves having subharmonic 
frequency and that oblique breakdown might be another limiting case with one secondary wave having 
zero frequency. Moreover, theoretical considerations suggested that oblique breakdown might be a stronger 
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nonlinear transition mechanism for two-dimensional boundary layers at supersonic speeds than any other 
resonance triad. 

A nonlinear transition mechanism, however, can only be considered to be an important physical mecha- 
nism if it can initiate the entire transition process of a laminar boundary layer to turbulence. This issue was 
addressed by the final part of the series of numerical studies, 6, 7,9 where the entire transition path of oblique 
breakdown for a flat-plate boundary layer at Mach 3 was simulated. The DNS in these references, however, 
had one shortcoming, that is, the flow was assumed to be symmetric in the spanwise direction to reduce 
computational costs. In order to ensure that this assumption is indeed valid, a new DNS has been conducted 
in which the artificial symmetry restriction was removed. As a first part, the final conference paper will 
show results from this simulation and additional longer time-average from the original symmetric simulation 
of Mayer et aV (CASE 3). Simulations of the complete transition process as in Mayer et al. 7 provide a 
valuable and extensive database that can be utilized for the validation of several different engineering tools 
across the various stages of transition. The second part of the final conference paper will therefore compare 
DNS results of the early and late nonlinear stages to data from NASA’s LASTRAC code. This code is based 
on the nonlinear parabolized stability equations. 

II. Simulation Setup for the DNS 

The simulation setup for the DNS in this abstract follows Mayer et al. 6,7 Supersonic flow at Mach 3 
over a flat plate is investigated numerically using spatial DNS. The physical conditions of the simulations 
match the Princeton wind tunnel conditions: 13 the unit Reynolds number formed with the free-stream 
velocity and free-stream viscosity at the inflow was Re = 2.181 x 10 e m~ 1 and the free-stream temperature 
is = 103. 6K. In the following, details of the computational setup are given and then information on 
governing equations, boundary conditions, disturbance generation and numerical method are repeated from 
Mayer et al. e ’ 7 for completeness and convenience of the reader. 

A. Computational Setup 

The main parameters for both DNS discussed in this abstract are summarized in table ??. The same 
nomenclature as in Mayer et al. 6,7 is employed. Hence, CASE 3 denotes the main case from the older 
references while CASE 7 is the newest simulation. The overall resolution is exactly the same for both cases. 
The grid is clustered in the streamwise direction using a fifth-order polynomial and in the wall-normal 
direction using a third-order polynomial. The computational grid in physical space consists of a total of 
about 212 million grid points for CASE 3 and 324 million grid points for CASE 7. 


Table 1. Main simulation parameters that differ between cases. 


Parameter 

Unit 

CASE 3 

CASE 7 

Domain size: 
XL 

M 

1.145 

1.050 

Vh 

M 

0.030 

0.030 

Grid size: 
n x 

h] 

2757 

2101 

n y 

h] 

301 

301 

K 

H 

128 

256 

n z 

h] 

255 

512 

n x x n y x n z 

h] 

211.6E6 

323.8E6 

Grid resolution (at outflow): 
points per A^ 1 ’ 1 ^ 

~ 440 

~ 440 

FFT’s: 

Symmetry in z? 

H 

yes 

no 


The inflow of the domain for both cases is located at Xq = 0.258m downstream of the leading edge 
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of the plate, whereas the outflow ranges from approximately 13.1 (CASE 7) to 14.5 (CASE 3) streamwise 
wavelengths X x of the oblique fundamental disturbance waves in the linear regime. The domain height is 
chosen as y* H = 0.030m ~ 5 boundary layer thicknesses S (laminar) at the outflow, such that even with the 
large increase in boundary layer thickness caused by the transition process no turbulent flow structures reach 
the free-stream boundary. 

Time-harmonic disturbances with a fundamental frequency of about /* = 6.36 kHz (F = 3 x 10 -5 ) are 
introduced through a blowing and suction slot located between x* = 0.394m and x\ = 0.452m (22 — Xi « A^). 
A discrete wave pair of instability waves with equal but opposite wave angle is excited for all cases. The 
spanwise wavenumber of /3* = 211.52m” 1 for the oblique wave pair is chosen such that the generated 
instability waves experience strong amplification as predicted by LST throughout the entire computational 
domain. This spanwise wavenumber determines also the domain width of all simulations, i.e. 2 ,*,, = X* = 
2t t//3* = 0.03m. 


B. Governing Equations 


The physical problem is governed by the conservation of mass, momentum and total energy with appropriate 
boundary and initial conditions. The fluid was considered to be a perfect gas with constant specific heat 
coefficients. The equations were cast in non-dimensional form using an arbitrary reference length (the plate 
length L* = 0.7239 m) and the free-stream values of the flow quantities at the inflow boundary: 


dp d 


+ nr (w) = 0 1 


dt dx 


d P Ui 1 d 1 1 x \ n 

~dT + ~dx~ ^ pUiUj + ° ijP ~ T v) =0 > 

dE 0 

~gf + ~^ Et+ P\ u i ~ u i T H + Qj) = 0 • 


dxj 


The total energy E t and the viscous stress are defined as: 


E t =p 


T 


UkUk 


k(k — 1 )Ma 


2 J 


Tin — 


p f dui 


duj 2 duk 

dxi 3 y dxk 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 


Re \dxj 

The pressure p is computed using the equation of state and the heat flux ^ is obtained from Fourier’s law: 

pT p dT 


P = 


Qi = 


nMa 2 ’ (k - 1 )Ma 2 RePr dxi ’ 

with k = 1.4 and Pr = 0.71. The viscosity is calculated using Sutherland’s law 


p(T) = T* 


1 + € 

T + ^r 


( 5 ) 


(6) 


with C = 110.4AT . 

C. Boundary Conditions 

At the inflow, the conservative quantities p , pui and E t , obtained from the similarity solution of a compressible 
flat-plate boundary layer, are specified. The outflow is treated with a buffer domain technique 14 to avoid 
reflections of disturbance waves. The buffer domain starts at 23 and ends slightly upstream of 2 l with a 
length of xl — 23 ~ 0.5Aa; for all cases. At the free-stream boundary, all total flow quantities are separated into 
base-flow and disturbance quantities. For the base-flow quantities, a homogeneous von Neumann condition is 
applied whereas for the disturbance quantities an exponential decay condition is employed that was derived 
for compressible flow using linear stability considerations . 15 In the lateral direction, periodicity is assumed. 
At the wall, the no-slip and no-penetration conditions are used except for the disturbance slot (see below). In 
addition, for the base flow, the wall temperature is set to the adiabatic wall temperature of the corresponding 
laminar flow, i.e. the initial condition, whereas temperature fluctuations are assumed to vanish. 
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D. Disturbance Generation 


The flow is forced over the disturbance slot by prescribing a time-harmonic function for the fundamental 
spanwise Fourier mode of the v- velocity. During the startup of the simulation, the forcing amplitude A(0, t) 
is ramped up in time over one disturbance period. The velocity distribution v p over the blowing and suction 
slot has the shape of a dipole and it is represented by a fifth-order polynomial that is smooth everywhere 
including at the end points: 


x p was defined as 


v ( x p , y = 0, 0, t) = A (0, t ) v p ( x p ) cos (-ut + 6 P (0)) . 


x 


p ~ 


2x - (X2 + Xi) 

X 2 - Xl 


-1 < Xp < 1 . 


( 7 ) 

( 8 ) 


E. Numerical Method 


The governing equations are integrated in time by employing a fourth-order Runge-Kutta scheme. The 
spatial derivatives are discretized using formally fourth-order split-finite differences in the streamwise and 
wall- normal directions. 16 The spanwise direction is assumed to be periodic and therefore transformed into 
spectral space using Fast Fourier transforms. Moreover, the spanwise discretization is pseudo-spectral, i.e. 
all nonlinear terms in the governing equations are computed in physical space and then transformed back 
into spectral space. Two options are implemented for the Fourier Transformations: (i) All flow variables (i.e. 
w-velocity, ^-velocity, etc.) are assumed to be symmetric to the centerline, except for the spanwise velocity 
w, which is antisymmetric. Symmetric quantities are then transformed into Fourier space using a Fourier 
cosine transformation and antisymmetric variables (w) are transformed by a Fourier sine transformation. 
Thus, only one-half spanwise wave length A- has to be computed for this configuration, (ii) No symmetry is 
assumed and therefore, all variables are transformed using a full Fourier transformation. This option requires 
the computation of the entire wave length X z in spanwise direction. 

The Fourier transformations are based on the VFFTPK library, which can be downloaded from netlib 
(http://www.netlib.org/vfftpack/). According to this library and its implementation in the Navier-Stokes 
code, a Fourier cosine transformation into spectral space and its back transformation into physical space are 
given by 


physical^spectral: 


spectral— ^physical: 


1 

2 (n z - 1) 


n z — 1 

</>§ + 2 ^ cos 

i=i 



K- 1 

4>0 + 2 0k COS 

k= 1 



(9a) 


(9b) 


for k = 0, ...,K — 1 and l = 0, ....,n z — 1, respectively. represent the Fourier amplitudes for mode k. 
Moreover, n z indicates the number of grid points used for resolving the spanwise direction in physical space 
over the interval [0, (n z — l)Ac] with 


Ar = 


A z 

2 (n z - 1) 


( 10 ) 


and K represents the number of modes in Fourier space (for the simulations n z = 2K — 1). 

The Fourier sine transformation to spectral space and its back transformation into physical space are as 
follows 


physical^spectral: 

& = f(<t>Y k 


spectral— >physical: 

« = {*)] 



n z — 1 

5Z <^ sin 
1=1 



K-l 

-2 J2 ^fc sin 


k= 1 



(lla) 


(lib) 
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for k = 0 , K — 1 and l = 0, n z — 1, respectively. 

In contrast to a symmetric simulation where only one-half of the spanwise wave length has to be calcu- 
lated, an asymmetric simulation requires the entire spanwise wave length as computational domain. Hence, 
for symmetric simulations n z represents the number of grid points in one-half wave length, whereas for 
asymmetric simulations, this number depicts the grid points in one full spanwise wave length. In this case, 
the grid spacing in spanwise direction is therefore obtained from 


A2- = 



(12) 


The full Fourier transformation for an asymmetric simulation is implemented according to 


physical— ^spectral: 

<t> o 

K 

spectral^physical: 


2 n 2 

1 

n z 

1 

n. 


Y 


1=0 
n z — 1 


v f 2nkl \ 

S * C0 T—) 


1=0 
n z — 1 


v — > . f 2irkl 

e * sm (— 

i=o x 


K - 1 


fo + Y 


k = 1 L 


~ c { 2nkl\ ~ . f 2irkl 


(13a) 

(13b) 

(13c) 


(13d) 


with k = 0, ..., K — 1 and l = 0, ...., n z — 1. As for the symmetric case, K denotes the number of Fourier 
modes. The entire storage space for the Fourier modes is however 2K — 1 since the cosine modes and the 
sine modes have to be stored separately. More details on the numerical method as well as validation cases 
can be found in the theses of Harris 16 and von Terzi. 17 


III. Results and Discussion 

This section discusses some preliminary results. In order to confirm that CASE 3 from Mayer et al. 7 
represents a valuable database for the validation of different engineering tools across the various stages 
of transition one remaining issue has to be clarified. In CASE 3, the spanwise direction was assumed to 
be symmetric since so far in the literature, oblique breakdown has always been initiated by two symmetric, 
oblique, first-mode type instability waves. It is, however, not straight forward to determine if this assumption 
is indeed justified. Hence, a final simulation where the symmetry restriction is removed is presented in 
section niEl Furthermore, this section also shows some results from CASE 3 that are averaged over a 
longer time interval. The final section (section III© provides a preliminary comparison of PSE results 
obtained from NASA’s LASTRAC code to CASE 3 for the early nonlinear transition regime of oblique 
breakdown. 

A. New DNS results 

The first two figures show results from the new simulation CASE 7. Since there is no symmetry condition in 
the spanwise direction, the computational domain for CASE 7 is doubled in this direction (see also Tabled 
and section II©. The influence of asymmetric modes on oblique breakdown initiated by two oblique waves 
with exactly the same amplitude and phase is limited since these modes are only generated by the round-off 
error of the calculation. In CASE 3 from Mayer et al . 7 the streamwise position of the final breakup into 
small-scale structures denotes the location where all modes with frequency unequal to integer multiples of 
the forcing frequency are strongly amplified. In CASE 7, a similar behavior can be observed. At exactly 
the same streamwise position (where the breakup into small-scale structures occurs) the asymmetric modes 
start to be amplified as illustrated by figure©. This figure shows the streamwise velocity disturbance of the 
first higher Fourier mode in the spanwise direction for the sine and cosine modes. In CASE 3, only cosine 
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* 


x , m 

Figure 1. Contours of streamwise velocity u obtained from CASE 7 for the first higher Fourier mode in the 
spanwise direction: (a) sine mode, (b) cosine mode. 


modes were calculated for the streamwise velocity because of symmetry and all sine modes were set to zero. 
Hence, the amplitude values of the sine mode in figure^i provide a measure for the magnitude of asymmetry 
in CASE 7. Since the sine mode in figure is more than 10 orders of magnitude smaller than the cosine 
mode, CASE 7 remains symmetric even after the breakup into small-scale structures. It is however clearly 
visible that the asymmetric modes are strongly amplified downstream of this position and will eventually 
reach high amplitude values in the turbulent region. 

The results from figure Q] corroborate that the symmetry assumption in CASE 3 is justified. Thus, for 
the rest of this section, we concentrate on CASE 3. The time average in Mayer et al. 6 7 is calculated over only 
12 forcing periods. If this interval length for the time average is sufficient, is depicted in figured 


0.0031 


0.0030 


0.0029 


0.0028 



t/T, 


forcing 


Figure 2. Skin- friction coefficient for CASE 3 as a function of interval length for time-averaging indicated by 
the number of forcing periods Tf orc i ng at three different streamwise positions. 


Figure 2 demonstrates the skin- friction coefficient as a function of the interval length for the time average 
at three different streamwise positions. For the first two positions, the skin-friction coefficient does not change 
significantly when the interval length for the time average is increased. At the last position (x* = 1.104m), 
however, a longer time average is required. A similar conclusion can be drawn from figure|21 which illustrates 
the streamwise distribution of selected mean-flow properties from CASE 3 for two different time averages. 
The curves with 12 forcing periods as time-average interval are very close to the curves with 25.75 forcing 
periods. 

As can be seen in figures |21 and 0 the increase in interval length for the time average does not strongly 
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Figure 3. Streamwise development of selected mean-flow properties from CASE 3 in comparison to different 
values published in the literature for turbulent supersonic flow 18 ’ 19 and theoretical models: 20 (a) skin-friction 
coefficient Cf, (b) Reynolds number based on momentum thickness 0. 


alter the mean values for CASE 3. However, the impact on fluctuation quantities, as for example the r.m.s. 
values, will be more pronounced. Hence, the final paper will discuss in more detail fluctuation quantities 
obtained from several different averaging intervals at different streamwise positions. Moreover, the final 
paper will also give a more detailed discussion of CASE 7. 

B. PSE Comparison 

For this abstract, only a preliminary comparison between PSE results and DNS data is provided. Figure 0] 
shows the streamwise development of wall-normal maximum in streamwise velocity disturbance u' for selected 
Fourier modes. Symbols represent PSE results and lines DNS data. For mode [1, 1], the agreement between 
PSE and DNS is excellent while for the other modes a discrepancy is visible. This discrepancy is due to the 
different disturbance generation in both methods. Note that the notation [h, k] is used to identify a particular 
wave according to its frequency h and its spanwise wavenumber k. h denotes multiples of the fundamental 
frequency /* = 6.36 kHz and k multiples of the smallest spanwise wavenumber (3* = 211.52 to - 1 . For the 
final paper, a detailed PSE study will be presented, in which the influence of the forcing method will be 
addressed. 


IV. Conclusions and Future Work 

Transition to turbulence via the oblique breakdown mechanism was investigated for a supersonic flat-plate 
boundary layer at Mach 3. Our previous studies of the same case focused on the detailed documentation 
of the different transition stages and demonstrated that oblique breakdown can lead to a fully developed 
turbulent boundary layer. In these studies, however, the flow was assumed to be symmetric in the spanwise 
direction. The verification of that assumption was an important part of the results presented in this abstract. 
To that end, a new simulation without the assumption of spanwise symmetry was performed. The final paper 
will provide a detailed discussion on this subject. Furthermore, new time-averaged data from CASE 3 of 
our previous study was shown. In the final paper, the impact of the time-averaging length on fluctuation 
statistics will be illustrated. Finally, the full paper will provide a detailed comparison of our DNS data to 
PSE results from NASA’s LASTRAC code. 
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Figure 4. Streamwise development of wall-normal maximum in streamwise velocity disturbance u' for selected 
Fourier modes: lines DNS, symbols PSE. 
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Oblique breakdown in a supersonic flat-plate boundary layer is investigated using Direct 
Numerical Simulations (DNS) and Parabolized Stability Equations (PSE). This paper con- 
stitutes an extension to our previous studies of the complete transition regime of oblique 
breakdown. In these studies, the flow was assumed to be symmetric in the spanwise direc- 
tion. A new DNS has been performed where the symmetry condition was removed. This 
simulation demonstrates that the “classical” oblique breakdown mechanism initialized by 
two symmetric instability waves with equal disturbance amplitudes loses its symmetry late 
in the turbulent stage for a low-noise environment. Hence, for the streamwise extent of the 
computational domain in our studies, the symmetry condition is justified. Furthermore, 
new data from a longer time average of the original symmetric simulation of oblique break- 
down (CASE 3) are discussed. These data verify that a converged time average is reached. 
The final part of the paper focuses on a comparison of PSE results obtained from NASA’s 
LASTRAC code to the DNS results. This comparison corroborates that the nonlinear PSE 
approach can successfully predict transition onset and that despite the large amplitude 
forcing used to introduce the oblique mode disturbances in the DNS, the latter constitutes 
a generic reference case for oblique breakdown at Mach 3 and, therefore, can be used to 
validate reduced order models for the full transition zone. 


N omenclature 


Latin 


Greek 


A 

Disturbance amplitude, 

Ot-i 

Streamw. amplification rate, 

c f 

Skin friction coefficient, 

Qt r , /? 

Streamw., spanw. wave number, 

Cp 

Specific heat at const, pressure 


Boundary layer thickness, 

E t 

Total energy, 

Sij 

Kronecker delta, 

f 

Frequency, 

At 

Interval length for time average, 

F 

Normalized frequency, 

Az 

Grid spacing in z, 

i 

vTT, 

e 

Phase in streamwise direction, 

k 

Thermal conductivity, 

0 

Momentum thickness, 

K 

Spanwise resolution (spectral), 

7 

Ratio of specific heats, 

L 

Reference length, 

A 

Wavelength, 

M 

Mach number, 

P 

Dynamic viscosity, 


Number of points in x , y, z, 

V 

Kinematic viscosity, 

n t 

Number of saved timesteps, 

p 

Density, 

P 

Pressure, 

T ij 

Stress tensor, 
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Pr 

Prandtl number, 

Qi 

Heat flux vector, 

Q 

Q-criterion, 

Re 

Reynolds number based on L, 

Ree 

Reynolds number based on 0, 

Rx 

Local Reynolds number, 

Ui 

Velocity vector, 

t 

Time, 

T 

Temperature, 

P forcing 

Forcing period, 

Xi 

Coordinate vector, 

x,y,z 

Streamw., wall-normal, spanw. directions, 

x 0 ,x L 

Location of inflow, outflow, 

Xl,X 2 

Start, end of disturbance hole, 

X3 

Start of buffer domain, 

VH 

Domain height, 

Zw 

Domain width, 


4> A flow quantity, 

to Angular frequency, 

Superscripts 

* Dimensional value, 

' Disturbance value, 

Time-averaged quantity, 
s Anti-symmetric (sine mode), 

c Symmetric (cosine mode), 

In Fourier space, 

Subscripts 

[h,k] Modes in [time, spanw. direction], 
oo Free-stream value, 

w Wall value. 


I. Introduction 


Boundary layer transition has important aerodynamic design implications on supersonic and hypersonic 
vehicles due to the strong increase in aerothermal loads. During the design process of the Space Shuttle 
Orbiter, boundary layer transition was recognized to be a significant aerothermodynamic challenge. 1 For 
accurate and successful prediction of transition onset, the transition process must be better understood in 
order to provide the future design community reliable physical prediction models. 2 These physical models 
need to incorporate the main transition stages, namely the receptivity regime, the initial linear disturbance 
development and the early nonlinear regime leading to the final breakdown. 

While the theoretical models governing the receptivity regime and the linear disturbance development 
are well established, 3-5 the nonlinear regime has received less attention. The knowledge of the relevant 
nonlinear mechanisms is, however, mandatory for the accurate determination of transition onset. Previous 
investigations 6-8 of the nonlinear transition regime for a two-dimensional, supersonic boundary layer discov- 
ered two main nonlinear mechanisms, oblique breakdown and asymmetric subharmonic resonance. Recently, 
a series of numerical studies 9-14 using direct numerical simulations (DNS) focused on several unresolved 
issues related to both mechanisms. In the first part of these studies, 9,10 the authors were able to identify 
oblique breakdown in the experiments by Kosinov and co- workers, 7, 15-1 ' who investigated asymmetric sub- 
harmonic resonance initiated by a wave train in a flat-plate boundary layer at Mach 2. The discovery of 
oblique breakdown in these experiments is of great importance since a detailed, experimental study of this 
mechanism has not been performed yet. 

The second part of these studies 13 addressed the question whether oblique breakdown or asymmetric 
subharmonic resonance is a stronger nonlinear breakdown mechanism in supersonic boundary layers. In 
order to answer this question, Mayer et al. 13 investigated the early nonlinear transition regime initiated 
by a broad disturbance spectrum on a cone at Mach 3.5. To excite a wide range of disturbance waves, a 
wave packet was generated by a pulse through a hole on the cone surface. The disturbance spectrum in 
the frequency-azimuthal mode number plane of the wave packet exhibited traces of oblique breakdown and 
new features that could be explained with new resonant triads that are composed of three boundary layer 
modes with three different disturbance frequencies. Furthermore, it was shown that asymmetric subharmonic 
resonance can be understood as a limiting case of these new resonance triads with both secondary waves 
having subharmonic frequency and that oblique breakdown might be another limiting case with one secondary 
wave having zero frequency. Moreover, theoretical considerations suggested that oblique breakdown might 
be a stronger nonlinear transition mechanism for two-dimensional boundary layers at supersonic speeds than 
any other resonance triad. A nonlinear transition mechanism, however, can only be considered to be an 
important physical mechanism if it can initiate the entire transition process of a laminar boundary layer to 
turbulence. This issue was addressed by the final part of the series of numerical studies, 11,12,14 where the 


2 of[18J 


American Institute of Aeronautics and Astronautics 


entire transition path of oblique breakdown for a flat-plate boundary layer at Mach 3 was simulated. Similar 
computations for first-mode type waves in a Mach 4.5 boundary layer were reported earlier by Jiang et al. ls 

The present paper complements the study by Mayer et a /. 11,12,1 4 via one additional DNS computation. 
In the earlier investigations, Mayer et aZ . 11, 12,14 assumed symmetry with respect to the spanwise direction 
to reduce computational costs. In order to ensure that this assumption is indeed valid, a new DNS has been 
conducted in which the artificial symmetry restriction was removed. Furthermore, the original symmetric 
simulation of Mayer et al . 11 ' 12,14 was continued and additional time-dependent data were saved in order to 
extend the interval for the time average and, hence, to ensure the convergence of statistical quantities. 

A major theme of this paper focuses on a detailed comparison between results obtained from a parab- 
olized stability equation (PSE) approach using the Langley Stability and Transition Analysis Codes (LAS- 
TRAC ) 19,20 and the DNS results. DNS results of the complete transition process as discussed in this paper, 
together with Mayer et al . 11,12 and von Terzi et al., 14 constitute a valuable and extensive database that can 
be utilized for the validation of several different engineering tools across the various stages of transition. With 
the development of LASTRAC, NASA Langley Research Center aims to create an integrated tool kit that 
can be employed for conventional N-factor calculations (based on LST or PSE) and more sophisticated sim- 
ulations from the receptivity process to the early nonlinear transition stages. Thus, a successful comparison 
between LASTRAC and the DNS results would validate LASTRAC’s nonlinear prediction capabilities. 


II. Governing Equations 


Supersonic flow at Mach 3 over a flat plate is investigated numerically using spatial DNS, linear, and 
nonlinear PSE. The governing equations are derived for a rectangular coordinate system with x as streamwise, 
y as wall-normal, and z as spanwise coordinates. The physical conditions of the simulations match the 
Princeton wind tunnel conditions . 21 The unit Reynolds number based on the free-stream velocity and 
free-stream viscosity is Re = 2.181 x 10 6 ?n _1 and the free-stream temperature is = 103. 6AT. Note 
that * indicates dimensional values. The fluid is considered to be a perfect gas with constant specific heat 
coefficients. The flow quantities are nondimensionalized by their approach-flow values, indicated by the 
subscript oo, except for the pressure and the total energy, which are scaled by the dynamic pressure pl^U^ 2 . 
The flow evolution is governed by the conservation of mass, momentum, and total energy: 

f + s; ( ^ > = 0, (1) 


d P Ui , d t , A- x n 

~dF + ~dxT 1 p * j + ° ijP ~ T v) = 0 > 
dE d 

~dT + dx ~^ Et+ p } Uj ~ UlTlj + ^ = 0 ’ 

where the symbols p and rq denote the fluid density and the velocity vector, respectively. 
The total energy E t and the viscous stress t,,- are defined as: 


E t = p 


/ T u k u k \ 


( dui 

duj 

2 r du k \ 

V7(7-l )M 2 2 ) 

1 ^ 
la; 

1 

^5 

\dxj 

dxt 

3 13 dx k ) 


( 2 ) 

(3) 

(4) 


with T as temperature. The pressure p is computed using the equation of state and the heat flux qi is 
obtained from Fourier’s law: 


pT p dT 

7 A I 2 ' Pl (7 — 1 )M 2 RePr dxt 


The viscosity is calculated using Sutherland’s law 


p(T) = rl 


1 + 


c 

T * 


T+^r 


( 6 ) 


with C = 110.4AT. Bulk viscosity is negligible. 
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The nondimensionalization of the governing equations introduces the Mach number M, the Prandtl 
number Pr, the ratio of specific heats, and the Reynolds number Re. The first three nondimensional 
parameters are defined as 


M = — -p- = —? U °° = 3, Pr = ^^= 0.71, and 7 = 1.4. (7) 

with c* poo , k* , and a*^ being the specific heat at constant pressure, the thermal conductivity, and the speed 
of sound of the approach flow, respectively. The Reynolds number, 


Re 


* nt r* 

P oc^oo-^ 


( 8 ) 


is based on a reference length L* , which is a constant for both numerical approaches. For the DNS, the plate 
length in the experiment is used while for the PSE, the similarity boundary-layer length scale ( L * = » / ) 


at the initial location Xg is chosen, 


L* DN s = 0.7239 and 


L 


PSE 


U* 


(9) 


A. Modifications for the PSE Approach 

Eqs. JIJ to ® are solved directly in the DNS. For the PSE approach, the governing equations are further 
simplified as discussed in detail by Chang . 19 All flow quantities </> are decomposed into the laminar mean 
flow solution $ and a disturbance fluctuation q V 


cf) = $ + <j)' . 

The disturbance fluctuation can be represented by the following wave ansatz 


</>' = e</> (x, y) exp 


a (£) + (3z — wt 


(10) 


( 11 ) 


with a = a r + iai as complex streamwise wave number, (3 as spanwise wave number, and u> as angular 
frequency. The integral form of the streamwise wave number is used to allow for a streamwise variation of 
a and to record the history effect. The shape function <j> is not only dependent on the wall-normal direction 
y as for linear stability theory (LST), but also on the streamwise coordinate x. Thus, non-parallel effects 
can be captured in the PSE approach through the variation of the streamwise wave number and the shape 
function. Substituting Eqs. itToll and ITTll into Eqs. GJ to © and subtracting the governing equations for 
the mean flow, leads to the final equation for the PSE approach 


L<) = /, with L = Ad x + Bd y + D — U yy dyy . 


(12) 


Here, the diffusion terms with respect to the streamwise direction x are neglected. The forcing function / 
contains the Fourier transform of the forcing applied in the PSE analysis and also the nonlinear terms for 
the nonlinear PSE approach. The coefficient matrices A, B , and C are functions of all wave numbers and 
the mean flow and can be found in Chang . 19 

For completeness, the frequency ui in Eq. CEH> is usually rescaled in stability calculations according to 


F = — 
R x ’ 


(13) 


where the Reynolds number R x is based on the similarity boundary-layer length scale L* 



R 


x — 



(14) 
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III. Simulation Setup for the DNS 


The simulation setup for the DNS follows Mayer et al . 11 ’ 12 and is repeated here for the convenience of the 
reader. In the next sections, the computational setup, the numerical method, boundary conditions and the 
disturbance generation is introduced. The purpose of the detailed description of the simulation setup, and 
especially the focus on the spanwise Fourier transformation in the numerical method section, is to provide 
the reader with all information necessary for the repetition of the simulations. Note that there is also a note 
on post-processing that discusses the temporal Fourier transformations applied to the DNS data to obtain 
the spectral composition of the disturbance held. 

A. Computational Setup 

The main parameters for both DNS discussed in this paper are summarized in Tabled The same nomen- 
clature as in Mayer et al. 11,12 is employed (see also nomenclature). Hence, CASE 3 denotes the main case 
from the older references while CASE 7 is the new asymmetric simulation. The streamwise and wall-normal 
resolution is exactly the same for both cases. The grid is clustered in the streamwise direction using a 
fifth-order polynomial and in the wall-normal direction using a third-order polynomial. The computational 
grid in physical space consists of a total of about 212 million grid points for CASE 3 and 324 million grid 
points for CASE 7. The difference in the overall number of grid points between both cases is mainly due to 
the slightly shorter domain in streamwise direction and the larger domain in spanwise direction for CASE 7. 


Table 1. Grid and simulation characteristics. Note that the notation [h,k\ is used to identify a particular 
wave according to its frequency h and its spanwise wave number k. h denotes multiples of the fundamental 
frequency and k represents multiples of the smallest spanwise wave number. 


Parameter 

Unit 

CASE 3 

CASE 7 

Domain size: 
XL 

[m] 

1.145 

1.050 

Vh 

M 

0.030 

0.030 

Grid size: 
n x 

h] 

2757 

2101 

n y 

h] 

301 

301 

K 

h] 

128 

256 

n z 

h] 

255 

512 

n x x n y x n z 

h] 

211.6E6 

323.8E6 

Grid resolution (at outflow): 
points per aL 1 ’ 1 ^ 

~ 440 

~ 440 

FFT’s: 

Symmetry in z? 

h] 

yes 

no 


The computational setup is illustrated in Fig. ^ The inflow of the domain for both cases is located at 
Xq ~ 0.258?n downstream of the leading edge of the plate, whereas the outflow ranges from approximately 
13.1 (CASE 7) to 14.5 (CASE 3) streamwise wavelengths A x of the oblique fundamental disturbance waves 
in the linear regime measured from the inflow. The domain height is chosen as y* H cs 0.030?n ss 5 boundary 
layer thicknesses 5 (laminar) at the outflow, such that even with the high increase in boundary layer thickness 
caused by the transition process no turbulent flow structures reach the free-stream boundary. 

Time-harmonic disturbances with a fundamental frequency of about /* = 6.36 kHz (F = 3 x 10 -5 ) are 
introduced through a blowing and suction slot located between x{ ~ 0.394m and x % — 0.452m (x 2 — x± ~ A^,). 
A pair of instability waves with equal but opposite wave angle is excited via the slot. The spanwise wave 
number of f3* ~ 211. 52m~ 1 for the oblique wave pair is chosen such that the generated instability waves 
experience strong amplification as predicted by LST throughout the entire computational domain. This 
spanwise wave number also determines the corresponding domain width for the simulations, i.e. z^y = A* = 
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2n//3* ~ 0.03m. For CASE 3, only half of this spanwise domain width is simulated because of symmetry, 
while, for CASE 7, the entire domain must be considered. 



B. Numerical Method 

The governing equations are integrated in time by employing a fourth-order Runge-Kutta scheme. The spa- 
tial derivatives are discretized using formally fourth-order split-finite differences in the streamwise and wall- 
normal directions. The spanwise direction is assumed to be periodic and, therefore, transformed into spectral 
space using Fast Fourier Transforms (FFT). Moreover, the spanwise discretization is pseudo-spectral, 22 i.e. 
all nonlinear terms in the governing equations are computed in physical space and then transformed back 
into spectral space. Two options are available for the spanwise Fourier transformations: (i) All flow variables 
(i.e. streamwise velocity u, wall-normal velocity v. etc.) are assumed to be symmetric to the centerplane, 
except for the spanwise velocity w, which is antisymmetric. Symmetric quantities are then transformed into 
Fourier space using a Fourier cosine transformation and antisymmetric variables ( w ) are transformed by 
a Fourier sine transformation. Thus, only one-half spanwise wave length X z has to be computed for this 
configuration, (ii) No symmetry is assumed and, therefore, all variables are transformed using a full Fourier 
transformation. This option requires the computation of the entire spanwise wave length X z . 

The Fourier transformations are based on the VFFTPK library, which can be downloaded from netlib 
(http://www.netlib.org/vfftpack/). According to this library and its implementation in the Navier-Stokes 
code, a Fourier cosine transformation into spectral space and its inverse transformation into physical space 
are given by 


physical— >spectral: 

4>l 

spectral^physical: 


1 

2 (n z - 1) 


n z — 1 

<t > o + 2 $ cos 

1=1 



K - 1 

$> + 2 ^ cos 

*;= i 



(15a) 


(15b) 


for k = 0, ...,K — 1 and l = 0, — 1, respectively. represents the Fourier amplitude for mode k. 

Moreover, n z indicates the number of grid points used for resolving the spanwise direction in physical space 
over the interval [0, (n z — l)Az] with 


Az = 


A z 

2 (n z - 1) 


(16) 


and K represents the number of modes in Fourier space (for the simulations n z = 2K — 1). 

The Fourier sine transformation to spectral space, with <fi s k as the Fourier amplitude, and its inverse into 
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physical space are 


physical— ^spectral: 

fa 

spectral^physical: 

fa 


1 


(n 2 - 1 ) 


n z — 1 

fa sin 


1=1 



K - 1 

-2 fa sin 

k=l 



(17a) 


(17b) 


for k = 0 , K — 1 and l = 0, n z — 1, respectively. 

In contrast to a symmetric simulation where only one-half of the spanwise wave length has to be calcu- 
lated, an asymmetric simulation requires the entire spanwise wave length as computational domain. Hence, 
for symmetric simulations n z represents the number of grid points in one-half wave length, whereas for 
asymmetric simulations, this number depicts the grid points in one full spanwise wave length. In this case, 
the grid spacing in the spanwise direction is therefore obtained from 


Az = 



(18) 


The full Fourier transformations for an asymmetric simulation are 


physical— ^spectral: 

fa 

fa 

fa 

spectral— ^physical: 

fa 


1 

2 n 

1 

n z 

1 

n z 


n z — 1 

- y fa> 
z 1=0 
n z — 1 

y fa cos 
1=0 
n z — 1 

y fa sin 
1=0 



K-l r 


*>+y 


k= 1 


facos 



+ fakSin 



(19a) 

(19b) 

(19c) 


(19d) 


with k = 0, ..., K — 1 and l = 0, ...., n z — 1. As for the symmetric case, K denotes the number of Fourier 
modes. The entire storage space scales by 2 K — 1 since the cosine modes and the sine modes have to be 
stored separately. More details on the numerical method as well as validation cases for the Navier-Stokes 
code can be found in the theses of Harris , 23 von Terzi , 24 and Mayer . 25 


C. Boundary Conditions 

At the inflow, the conservative quantities p, piq, and E t , obtained from the similarity solution of a com- 
pressible flat-plate boundary layer, are specified. The outflow is treated with a buffer domain technique 26 
to avoid reflections of disturbance waves. The buffer domain starts at X 3 and ends at xl with a length of 
xl — x 3 ~ 0.5A X (Fig.Q. At the free-stream boundary, all total flow quantities are separated into base flow 
and disturbance quantities. For the base-flow quantities, a homogeneous von Neumann condition is applied 
whereas for the disturbance quantities an exponential decay condition is employed that was derived for com- 
pressible flow using linear stability considerations. 2 ' Periodicity is assumed in the spanwise direction. At 
the wall, the no-slip and no-penetration conditions are used except for the disturbance slot (see below). In 
addition, for the base flow, the wall temperature is set to the adiabatic wall temperature of the corresponding 
laminar flow, i.e. the initial condition, whereas temperature fluctuations are assumed to vanish. 


D. Disturbance Generation 

The flow is forced through the disturbance slot by prescribing a time-harmonic function for the fundamental 
spanwise Fourier mode of the ^-velocity. During the start-up of the simulation, the forcing amplitude A{0) 
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is ramped up in time over one disturbance period. The velocity distribution v p over the blowing and suction 
slot has the shape of a dipole and it is represented by a fifth-order polynomial that is smooth everywhere 
including at the end points 


v (x p , y = 0,/3,t) = A ((3) v p (x p ) cos (-wt + 8 P (/?)) . 


x p is defined as 


and v p as 


x p (x) = 


2x - (X 2 + Xi) 
x 2 - X\ 


for — 1 < x p < 1 , 


i{Xp) — 


L.5 4 (1 + x p ) 3 (3 (1 + x p ) 2 - 7 (1 + x p ) + 4) , x p < 0 


-1.5 4 (1 - x p f (3 (1 - x p ) 2 - 7 (1 - x p ) + 4 ) , x p > 0 . 

The amplitude A = 0.003 and, without any loss of generality, 6 p = 0 for the DNS described herein. 


( 20 ) 

( 21 ) 

(22) 


E. Post-processing 

In order to obtain the spectral composition of the disturbance field from the DNS in a post-processing step, 
the physical time signal is Fourier transformed in time. The post-processing tool is based on the EAS3 tool kit 
from the Universitat Stuttgart (http://en.wikipedia.org/wiki/EAS3). In EAS3, the Fourier transformations 
are given as follows 


physical— ^spectral: 

Po 


Pm 


1 

n t 


nt — 1 

'y ' Pn J 
n — 0 


1 

n t 


nt — 1 

y 2p n exp 

n — 0 


—i2nmn\ 
nt ) 


(23a) 

(23b) 


for m = 1, ..., M — 1 and M = nt, /2. Here, M contains the number of Fourier modes and nt the number of 
sample points of the physical signal. The output data of the tool is the absolute value \p m \ and the negative 
phase — arg (p m ) of the complex number p m . The back transformation is defined as 


spectral— ^physical: 


rM-i 


Re 


y, pm exp f 

kv) — n \ 


_m = 0 


( i2Trmn\ 
nt ) 


for n = 0, ..., nt — 1. 


(24a) 


IV. Simulation Setup for the PSE Approach 

In this section, the simulation setup of the PSE calculation is given. The computational domain differs 
from the DNS due to the different theoretical and numerical approaches. A short discussion on the numerical 
method, the boundary conditions and initial condition is also provided. 

A. Computational Setup 

The Reynolds number range and the disturbance frequency of the PSE calculations and of the DNS (CASE 3) 
are depicted in the stability diagram of Fig. |3 Two main PSE calculations have been performed, which 
vary in the starting position for the marching procedure. For one PSE calculation, from hereon abbreviated 
as PSE 1, the starting position is located close to the first neutral branch at Xq SE1 = 0.068 m ( R x ~ 385) 
whereas, for the second PSE calculation, PSE 2, the domain starts at x ESE2 = 0.6 m ( R x ~ 1144). The 
reasoning behind these choices is explained in section V|2j The domain height is chosen as 200 x L *, with L* 
being the similarity boundary-layer length scale. The boundary layer is disturbed by prescribing non-parallel 
eigenfunctions obtained from a precursor non-parallel stability analysis (for a frequency of f* = 6.36kHz 
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Figure 2. Stability diagram for f3* ~ 211.52m 1 used to illustrate the computational domain for the PSE 
calculations. The box between R x = 927 and R x = 993 indicates the position of the blowing and suction slot in 
the DNS. 


and a spanwise wave number of /3* ~ 211.52m -1 ) at the inflow. The spanwise direction is assumed to be 
symmetric and is discretized using Fourier transformations. The number of Fourier modes for the oblique 
breakdown computations is selected to ensure an O(10 -4 ) decay in amplitude from the most energetic modes 
to the tail of the disturbance spectrum. The wall-normal direction is resolved by 171 points with clustering 
of grid points near the surface as well as the critical layer of the first-mode type waves. 

B. Numerical Method 

The governing linear and nonlinear PSEs are solved by a fourth-order central difference scheme in the wall- 
normal direction and by a first- or second-order one-sided difference scheme in the streamwise direction. The 
spanwise direction is assumed to be periodic and therefore transformed into spectral space using Fourier 
transforms. Furthermore, the spanwise discretization is pseudo-spectral, 22 i.e. all nonlinear terms on the 
right-hand side of Eq. C0 are computed in physical space and then transformed back into spectral space. As 
for the DNS, the spanwise direction can be set to be symmetric in it-velocity, v- velocity, temperature, density, 
and pressure and antisymmetric in wi-velocity reducing the computational domain to half a fundamental 
spanwise wave length A*. This approach is adopted for the present paper. Details on the discretization 
scheme can be found in Chang. 28 

C. Initialization of PSE Calculations and Boundary Conditions 

To initialize the PSE analysis, mean-flow data for all streamwise positions and the disturbance information 
at the starting location of the streamwise marching procedure are required. A convenient way to initialize 
the PSE analysis is to prescribe eigenfunctions from LST. This approach, however, would result in a “tran- 
sient” region where the PSE solution transitions from the LST eigenfunctions, satisfying parallel theory, at 
the inflow to the solution of the non-parallel parabolized stability equations. Hence, in the present work, 
eigenfunctions obtained from a precursor non-parallel stability analysis are used. The mean-flow data stems 
from a compressible similarity solution with adiabatic wall boundary condition. The differences between the 
similarity solution and the DNS mean flow were found to be negligible. At the wall, the no-slip condition 
is enforced and temperature fluctuations are set to zero while at the free stream, non-reflective boundary 
conditions are employed. 
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V. Results and Discussion 


The results section is split into three main parts. The first part gives a summary of the main findings 
for CASE 3 in the transitional regime, which are similar to CASE 2 in Mayer et al. 11 Note that CASE 3 is 
symmetric in spanwise direction since, so far in the literature, oblique breakdown has always been initiated 
by two symmetric, oblique, first-mode type instability waves. Furthermore, the results of CASE 3 are 
complemented by discussing CASE 7, for which the symmetry restriction is removed. The second part of 
this section (section VO is focused on longer time averages for the turbulent regime and their comparison 
with shorter time averages from Mayer et al. 12 The final part of the results section (section V0 provides a 
comparison of PSE results obtained from NASA’s LASTRAC code to results from CASE 3 for the receptivity 
regime and the early nonlinear transition regime of oblique breakdown. This comparison provides additional 
validation of the DNS setup and the nonlinear PSE approach of LASTRAC. 

A. Removing the Symmetry Condition in the Spanwise Direction 

Before the influence of the symmetry condition on CASE 3 is evaluated using results from the new asymmetric 
simulation (CASE 7), some main findings for CASE 3 from the transitional regime are summarized here. As 
already mentioned in Mayer et al. 11 and in von Terzi et al., 14 the initial growth of the forced modes [1, ±1] 
starts to deviate from linear growth when the higher-harmonic inodes, as for example the stationary modes 
[0, ±2] or modes [1, ±3], reach comparable amplitude levels to modes [1, ±1] (see Fig. |3Ji). This event marks 
the end of the early nonlinear regime (at about x* = 0.7 to). Note that the notation [h, k] is used to identify 
a particular wave according to its frequency h and its spanwise wave number k. h denotes multiples of the 
fundamental frequency and k represents multiples of the smallest spanwise wave number. 




Figure 3. Streamwise amplitude development of mode shape maximum for modes with fundamental frequency 
(a) and subharmonic frequency from CASE 3 (b): For (a): (o, — ) [l,d=l], (□, — ) [1,=L3], (0? -•-) [l,=b5], (A, 
[l,d=7]); For (b): ( — ) [0.5,0], (- -) [0.5, ±1], [0.5, ±2], [0.5, ±3]. Results in (a) are from CASE 1 and 2 in 

Mayer et al. 11 and in (b) from CASE 3. 

At about x* = 0.8 m, nonlinear saturation of the u - velocity disturbance amplitude sets in (Fig. [3). Modes 
that are not directly generated by nonlinear wave interactions from the initial forced modes [1,±1], as for 
example modes with subharmonic frequency in Fig. 0), also experience amplitude growth and saturation at 
about x* = 0.8 m. These modes are initialized by round-off errors due to the limited machine precision. One 
major difference between Fig. IHK and Fig. 0) is the sudden increase in the streamwise growth rate for all 
modes with subharmonic frequency at about x* = 0.9 m. This sudden increase in amplification seems to be 
linked to the breakup of large longitudinal structures into small-scales as depicted in Fig. Q] (between about 
x* = 0.84m and x* = 0.9m). This figure shows instantaneous three-dimensional isosurfaces obtained from 
the Q-criterion 29 (Q=15000). Q is related to the second invariant of the velocity gradient tensor. Positive 
values of Q reveal flow regions where rotation dominates the strain. 

The flow structures in Fig. are symmetric with respect to the centerplane of the domain. This is to 
be expected since these figures are obtained from CASE 3, where symmetry is enforced. If the symmetry 
condition is not enforced as for CASE 7, the picture does not change as demonstrated by Fig. [3 This figure 
illustrates contours of streamwise velocity u of the first higher Fourier mode (k = 1) in spanwise direction 
from CASE 7 for the sine and cosine modes (Eq. 1 1911 . respectively. The minimum and maximum of the 
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(a) 



Figure 4. Instantaneous flow structures identified by the Q-criterion for Q = 15000 (CASE 3) between 
x* = 0.798m and x* = 0.924m. Also shown are contours of spanwise vorticity at z* ~ —0.0087m. (a) Entire 

three-dimensional view, (b) close-up of the breakdown region confirming that the Q-criterion predicts similar 
structures as illustrated by the spanwise vorticity; M=3.0, TJ o =103.6K, flat plate. 


contour levels in Figs. EH and b are different in order to emphasize the flow structures. The influence of 
asymmetric modes on oblique breakdown initiated by two oblique waves with exactly the same amplitude 
and phase is limited since these modes are only generated by the round-off error of the calculation. For 
CASE 3, the streamwise position of the final breakup into small-scale structures denotes the location where 
all modes with frequency unequal to integer multiples of the forcing frequency are strongly amplified (at 
about x* = 0.9m in Fig. EH). For CASE 7, a similar behavior can be observed. At the same streamwise 
position (where the breakup into small-scale structures occurs) the asymmetric modes also start to become 
amplified as illustrated by Fig. EH. The amplitude values of the streamwise velocity for the sine mode in 
Fig-13 1 provide a measure of the magnitude of asymmetry in CASE 7. For the it- velocity, this mode is set to 
zero in CASE 3. Since the contour levels for the sine mode in Fig. Eli are more than 10 orders of magnitude 
smaller than the contour levels for the cosine mode, CASE 7 remains essentially symmetric even after the 
breakup into small-scale structures. This is true over the entire domain length of CASE 3. It is however 
evident that the asymmetric modes are strongly amplified in the downstream direction and will eventually 
reach high amplitude levels. 

B. Longer Time- Averaging 

The early turbulent regime for oblique breakdown (CASE 3) is discussed in detail in Mayer et al. 12 It is 
shown how the skin friction coefficient increased in the streamwise direction during transition until a peak 
was reached. Downstream of the peak, a loss of periodicity in the time signals for the investigated setup 
occurred and the skin-friction approached correlations and comparable data for turbulent boundary layers 
in the literature. Furthermore, a logarithmic region in the van Driest transformed mean streamwise velocity 
profile was detected and power spectra of the velocity components matched well-known theoretical scaling 
laws. 

In Mayer et al., 12 turbulent quantities are calculated from data, which are time- averaged according to 

ii r^z pto+At 

<(>= — At/ J <t> (t, z) dtdz . (25) 

Moreover, fluctuations about the mean of a Reynolds- averaged quantity are denoted by 4>' and fluctuations 
about the mean of a Favre-averaged quantity are denoted by 4>" . The interval length At is, however, 
only 12 forcing periods Tf orc i ng in Mayer et al. 12 Whether this interval length is sufficiently long, is discussed 
in this section. Fig. El demonstrates the skin-friction coefficient as a function of the interval length for the 
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Figure 5. Contours of instantaneous streamwise velocity u obtained from CASE 7 for the first higher Fourier 
mode in spanwise direction: (a) sine mode, contour levels from — 1.0E— 12 to 1.0E— 12, (b) cosine mode, contour 
levels from —0.1 to 0.1. 
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Figure 6. Skin-friction coefficient for CASE 3 as a function of interval length for time-averaging indicated by 
the number of forcing periods Tf orc i ng at three different streamwise positions. 


time average at three different streamwise positions. These positions are downstream of the skin-friction 
peak, within the early turbulent region, where the flow starts to lose its temporal periodicity with respect to 
the forcing frequency. For the first two positions, the skin-friction coefficient does not change significantly 
when the interval length for the time average is increased. This is most likely due to the fact that the flow 
field still exhibits some temporal periodicity and a negligible low-frequency disturbance content. Hence, 
a time interval of 12 forcing periods seems to be sufficient for a good time average. At the last position 
(a;* = 1.104m), however, a longer time average is required (even longer than A t/Tf orcing = 25.75). Time 
signals from this position are shown in Mayer et al. 12 and they demonstrate that the flow field has lost its 
temporal periodicity at this location. A similar conclusion can be drawn from Fig.0 which illustrates the 
streamwise distribution of selected mean-flow properties from CASE 3 for two different time averages. The 
curves with 12 forcing periods as time-average interval are very close to the curves with 25.75 forcing periods. 

As can be seen in Figs. ED and 0 the increase in interval length for the time average does not strongly 
alter the mean values for CASE 3. However, the impact on fluctuation quantities, as for example the r.m.s. 
values, might be more pronounced. Therefore, Fig. 0 compares the r.m.s. values of the streamwise velocity 
component and the Reynolds shear stress u"v" obtained from data with At/Tf orc i ng = 12 to the data with 
A t/Tf orc i n g = 25.75 at x* = 1.087m. In addition, profiles published in the literature from the temporal DNS 
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Figure 7. Streamwise development of selected mean-flow properties from CASE 3 in comparison to different 
values published in the literature for turbulent supersonic flow 30 ’ 31 and theoretical models: 32 (a) skin-friction 
coefficient Cf , (b) Reynolds number based on momentum thickness ©. 


by Guarini et al. 30 of a compressible turbulent boundary layer at Mach 2.5 are also included in this figure. 
The profiles from the different time averages are in surprisingly good agreement leading to the conclusion 
that even for the fluctuation quantities, a time interval of A t/Tf orc i ng = 12 for the average is satisfactory 
up to x* = 1.087 m. 



Figure 8. Wall-normal distribution of r.m.s. values in streamwise direction (a) and Reynolds shear stress u"v" 
(b) at x* = 1.087m for CASE 3. 


C. PSE Comparison 

As a short summary, the first two sections (section V^l and V[Bl confirmed that CASE 3 from Mayer 
et al. 12 is a representative case for oblique breakdown at Mach 3. The symmetry condition is applicable for 
the computational setup of CASE 3 and the data was sufficiently time-averaged. In this section, we employ 
NASA’s LASTRAC code in order to study its nonlinear PSE capabilities and to evaluate the simulation 
setup of CASE 3 in the receptivity regime and early nonlinear regime. 

1. Receptivity Study 

In order to investigate the receptivity behavior of a laminar Mach 3 boundary layer with adiabatic wall, 
LASTRAC’s receptivity module is applied. This module can be used for making integrated predictions of 
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receptivity and subsequent evolution of instability waves for certain simple classes of receptivity mechanisms, 
including the receptivity to unsteady blowing and suction at the surface (as in CASE 3). Before presenting 
any results for the instability wave perturbations excited by the unsteady blowing and suction slot from 
Eq. (®, we intend to explain the receptivity behavior of a supersonic boundary layer to this form of forcing 
in general. As described in the context of the DNS results, a localized blowing and suction slot excites 
the first mode-type instability wave at the frequency (and spanwise wave number) of excitation as well as 
other disturbances that are stable and, hence, do not influence the disturbance dynamics sufficiently far 
downstream. We confine our attention to the instability wave portion of the boundary layer response and 
examine how well the theoretical models (and, in particular, the model based on adjoint PSE 20 ) can predict 
the first-mode type perturbations in this specific case. The theoretical prediction for the instability wave 
perturbations due to the surface actuator can be expressed in the form of Eq. <HU), where 

e = -4= f A (£, p) F (£, (3) exp (-*6> r (£, /?)) . (26) 

V "7T J xq 

Here, the efficiency function A represents the instability wave portion of the Green’s function associated with 
the boundary layer response to this type of wall forcing, i.e. , the effective initial amplitude of the instability 
wave due to a point source actuation. It characterizes the intrinsic receptivity characteristics of the local 
boundary layer to a specific type of forcing and is independent of the spatial distribution of the source. The 
dependence on the source geometry is reflected via the convolution of the efficiency function A with the 
spanwise Fourier transfrom of the unsteady blowing and suction velocity v p from Eq. in 

F(x,(3) = ^Av p (x) . (27) 

The quantity 0 r in Eq. m denotes the phase of the velocity perturbation 

6r(x,/3)=[ a r (£,/3)d£. (28) 

J Xq 




Figure 9. Receptivity behavior of a laminar Mach 3 boundary layer with adiabatic wall due to blowing and 
suction according to Eq. applied in CASE 3: (a) efficiency function A from PSE as a function of streamwise 
position x, and (b) predicted amplitude value for modes [1, =bl] from PSE (solid line) and from DNS (dot) at 
position x* = 0.5 m as a function of starting position of disturbance slot x\. Note that the inset in (b) shows 
the combined amplitudes of modes [=Ll,d=l]. 


Fig.|2lillustrates the magnitude of the efficiency function A, plotted as a function of the streamwise source 
location x. The choice of normalization of the efficiency function is the same as in Chang and Clioudhari. 20 
It can be seen in Fig. that the intrinsic efficiency values vary weakly with the source location, with less 
than 33% variation over a broad range of actuator locations. Thus, the dominant effect of source location 
on the amplitude of modes [1, ±1] is due to the progressive reduction in the linear amplification potential 
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associated with source locations downstream of the first neutral branch. This effect is reflected in the strong 
reduction in the amplitude of modes [1, ±1] at x* = 0.5m as the starting location for the blowing and suction 
slot is varied along the streamwise direction in Fig. |2l)- Note that, for the results plotted in this figure, the 
measurement location is sufficiently downstream of the source, but still within a predominantly linear region 
in terms of the fundamental mode evolution. The amplitude of the instability wave extracted from the 
DNS computation (shown as a dot in Fig. |5J}) is within 3% of the value predicted by the theory. This 
small discrepancy could be even further narrowed by accounting for the cumulative effects of disturbance 
nonlinearity between the disturbance slot and the measurement location, but was deemed unnecessary. 

2. Finite Amplitude Forcing 

The early nonlinear regime of oblique breakdown at Mach 3 is studied by performing PSE calculations with 
two different simulation setups. As already mentioned in section ivm the first setup starts at the first 
neutral branch (PSE 1) while the second PSE calculation starts at x* = 0.6 m (PSE 2). For PSE 1, the 
initial disturbance amplitudes of modes [1, ±1] at the neutral branch are scaled in order to match the absolute 
amplitude values from the DNS farther downstream for this mode, while for PSE2, the initial amplitudes of 
the first 6 dominant modes ([1, ±1], [0, ±2], and [1, ±3]) are scaled to the corresponding values from the DNS 
at x* = 0.6m. The reasoning for the setup PSE 1 is to examine whether the DNS results (CASE 3) can be 





Figure 10. Streamwise amplitude development of the wall-normal maximum (mode shape maximum) for 
selected modes calculated from physical data according to the Fourier transforms given in Eqs. ltl5ll and (f25b : 
(a) (symbols) PSE 1, (dashed lines) DNS, (b) (symbols) PSE 2, (dashed lines) DNS. (c) Setup from PSE 1, 
but with (solid lines) and without (symbols) the suppression of modes [0, d=2]. 
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completely recovered by a PSE calculation starting with a very low disturbance amplitude at the first neutral 
branch, which would be the case during natural transition (i.e. without an artificial disturbance source at 
the plate surface as employed in the DNS). As shown in Fig. II Oh,, which illustrates the streamwise amplitude 
development of the wall-normal maximum (mode shape maximum) for selected modes, the DNS cannot be 
completely recovered resulting in a discrepancy between the nonlinear generated modes from the DNS and 
the nonlinear PSE analysis. This difference is mainly caused by nonlinear effects during the receptivity 
process in the DNS since the forcing amplitude of the DNS is already quite large (A=0.3%). If in a second 
PSE analysis the setup is changed to PSE 2 so that the amplitudes of not just the fundamental, but also 
the most dominant nonlinearly excited modes are matched with the DNS, then the agreement between the 
data from PSE and the data from DNS improves significantly as demonstrated in Fig. 110b. 

We also examined the effect of suppressing modes [0, ±2] to disrupt the energy cascade from modes 
[1,±1] to their higher-harmonic modes associated with oblique breakdown. Mode [0, ±2] was canceled out 
by a modest stationary forcing at x* = 0.6m in a new PSE calculation for setup PSE 1. Fig. UU' shows 
the results corresponding to this simulation. Symbols refer to the PSE results from Fig. II Oh ., where modes 
[0, ±2] are not suppressed and solid lines represent the new simulation. Clearly, the influence of nonlinear 
wave interactions is weakend resulting in a delay of transition. While the delay in transition onset location 
is approximately 5% increased, it may be less in a natural disturbance environment. 
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Figure 11. Streamwise development of skin friction coefficient: (a) (solid line) DNS, (symbols) PSE, (b) PSE 
with (solid line) and without (dashed line) suppression of modes [0,=b2]. 

The previously mentioned findings are confirmed by the streamwise development of the skin friction 
coefficient in Fig. EH The comparison between DNS and PSE results (PSE 2) in Fig. II Ih leads to an 
excellent agreement between both methods indicating that nonlinear PSE is able to predict the transition 
onset correctly. In fact, the agreement extends up to the location, where the skin friction coefficient is nearly 
twice the value of the underlying laminar skin friction. This suggests that (i) the PSE results could, perhaps, 
be bridged with a turbulence model to provide integrated predictions of transition and turbulence 33 and that 
such PSE results could be used, as a cost effective means, to initiate simulations of fully turbulent boundary 
layers. For the case with suppressed modes [0,±2], the transition onset is indeed moved downstream as 
depicted in Fig.HTI). 


VI. Conclusions 

Transition to turbulence via oblique breakdown in a Mach 3 boundary layer was investigated using 
DNS and PSE. Our previous studies of the same case focused on the detailed documentation of the different 
transition stages and demonstrated that oblique breakdown can lead to a fully developed turbulent boundary 
layer. In these studies, however, the flow was assumed to be symmetric in the spanwise direction. A 
new DNS was performed where the symmetry condition was removed. This simulation demonstrated that 
oblique breakdown, initialized by two oblique instability waves with exactly the same amplitude level, loses 
its symmetry late in the turbulent stage for a low- noise environment. Hence, for the streamwise extent of 
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the computational domain in this study, the symmetry condition is justified. Furthermore, new data from 
a longer time average of the original symmetric simulation of oblique breakdown (CASE 3) was discussed. 
These data confirm that a converged time average is reached. The final part of the paper focused on a 
comparison of PSE results obtained from NASA’s LASTRAC code to the DNS results all the way from 
the disturbance generation to the initial part of the breakdown process. Good agreement was found for all 
different transition stages, including receptivity, linear, and nonlinear disturbance evolution, and the initial 
part of the breakdown process until the mean wall shear had nearly doubled with respect to its laminar 
value. Thus, the present paper corroborates that our oblique breakdown case for Mach 3 constitutes a 
generic reference case that can be successfully used to validate different engineering tools, as for example 
the PSE approach, across the different transition stages. Furthermore, the PSE approach combined with 
turbulence models for transitional flow might be a sophisticated solution for modeling the entire transition 
process of a supersonic boundary layer. 
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